Protecting coherence in Optimal Control Theory: state dependent constraint 

approach 
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Optimal control theory is developed for the task of obtaining a primary objective in a subspace 
of the Hilbert space while avoiding other subspaces of the Hilbert space. The primary objective 
can be a state-to-state transition or a unitary transformation. A new optimization functional is 
introduced which leads to monotonic convergence of the algorithm. This approach becomes necessary 
for molecular systems subject to processes implying loss of coherence such as predissociation or 
ionization. In these subspaces controllability is hampered or even completely lost. Avoiding the 
lossy channels is achieved via a functional constraint which depends on the state of the system at 
each instant in time. We outline the resulting new algorithm, discuss its convergence properties 
and demonstrate its functionality for the example of a state-to-state transition and of a unitary 
transformation for a model of cold Rb2. 

PACS numbers: 



I. INTRODUCTION 

Coherent control utilizes the wave properties of mat- 
ter to steer a quantum dynamical process to a desired 
outcome. The source of control is interference, construc- 
tive to achieve the goal and destructive to eliminate un- 
wanted consequences [T], Q. The agents of control are 
external fields, in particular electromagnetic fields. The 
experimental and theoretical challenge lies in identifying 
these control fields. The present study is aimed at finding 
control fields which are constrained to limit the damage 
which the control field may cause to the controlled sys- 
tem. 

Theoretically, the control problem can be formulated 
as an inversion: finding the field subject to the quantum 
dynamics which leads to the desired outcome. Optimal 
control theory (OCT) has been developed as a tool to 
address this problem [H, H| . It can be formulated starting 
from a variational ansatz Q or using Krotov's method 
[H, 0] . Recently the Krotov's method has been extended 
to include a strict limitation on the spectrum of the op- 
timized field 0]- 

The most well studied task in OCT has been the goal 
of a state-to-state transition. Given an initial state ipi n i 
and a closed quantum system, the field needs to be found 
which drives the system to a specific final state ipfi n - 
This task has been shown to be completely controllable 
HQ if the fields are not restricted. Moreover the control 
landscape is favorable composed of flat ridges such that 
the climb in the gradient direction will lead to one of the 
many possible solutions [1(3] ■ 

A more involved control task is to optimize the ex- 
pectation value of an operator at a final time, (B(tjj n )). 
This task can be formulated in the framework of open 
quantum systems. The OCT approach yields an itera- 
tive solution to the inversion problem which is based on 
propagating the system density operator pg(t) forward 



in time and the target operator B(i) backward in time 
iil- 

The prospect of quantum computing has posed an even 
more complex control problem: imposing a unitary trans- 
formation U on a subset of quantum states which act as 
the quantum register. The unitary transformation car- 
ries out a specific computational task. This control task 
is equivalent to N simultaneous state-to-state transfor- 
mations @, [3 ■ The solution of the iterative set of 
equations has been shown to become exponentially more 
difficult with the size N of the unitary transformation 
Q. These findings are in accordance with a very com- 
plex control landscape [l6l |. 

A further step up in complexity is the task of imposing 
a unitary transformation under dissipative conditions. 
This task emerges in the quantum governor [T7| , in quan- 
tum information processing and it is a traditional task in 
nuclear magnetic resonance (NMR) spectroscopy [Tsl | . 

In any practical implementation the positive task of 
obtaining the final goal has to be weighted by possible 
negative consequences. For example control fields of high 
intensity can damage the system by causing ionization 
or dissociation. A remedy for this problem consists in 
restricting the population in certain lossy excited state 
manifolds. This task has been the motivation for the de- 
velopment of local control theory (LCT) [lj, Hj. LCT 
has been applied to lock unwanted electronic excitations 
[2ll | and recently to the problem of quantum informa- 
tion proces sing where avoiding population loss becomes 
crucial [22l. |23||. 

However, OCT is more powerful than LCT and it is 
therefore desirable to incorporate constraints describing 
negative consequences of the control process into the al- 
gorithm. Such constraints depend on the state of the sys- 
tem at intermediate times [2J, |25|, |26|, l27| . For example, 
the system can simply be restricted to remain in an "al- 
lowed" or to avoid a "forbidden" subspace during its evo- 
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lution. More elaborate examples include imposing: a pre- 
defined path between an initial and a final state [261 ] , and 
maximizing the expectation value of a time-dependent 
operator throughout the optimization time interval (25j . 

Previous OCT studies which impose state-dependent 
constraints [13, US were performed for state-to-state 
optimization and are based on the variational approach. 
In the present work, an optimization algorithm including 
state-dependent constraints is obtained using the Krotov 
method for the state-to-state case as well as for unitary 
transformations. The Krotov method offers the advan- 
tage that the monotonic convergence of the algorithm 
can be ensured by the choice of the imposed constraints 
0, HI]- A brief comparison with the algorithms using 
the variational ansatz [U HE, HH| will be given for the 
state-to-state optimization. 

The paper is organized as follows: The state-dependent 
constraints are formulated in Sec. [TH and the resulting 
algorithm is presented for optimization of state-to-state 
transition and of a unitary transformation. A review of 
the Krotov method together with an outline of the deriva- 
tion of the equations presented in Sec. [IT] is given in the 
Appendix. Sec. IIIII introduces a model example and il- 
lustrates optimization under state-dependent constraints 
for a state-to-state transition and for a unitary transfor- 
mation. Our findings are compared to related approaches 
in Sec.|TV] Finally, Sec. IVl concludes. 

II. FORMULATION OF STATE-DEPENDENT 
CONSTRAINTS IN THE KROTOV METHOD 

A. Optimization of a state-to-state transition 

The dynamics of the system is governed by the time 
dependent Schrodinger equation 

j t \ V {t)) = -ifi[e(t)M*)>, (1) 

where \<p(t)) represents the state of the system at time t, 
and 

H[e] = H - £e(i), (2) 

is the system+control Hamiltonian. Ho denotes the field- 
free Hamiltonian, e(t) the semiclassical control field and 
p, is a system operator describing the coupling between 
system and field. 

The objective of the optimization is to find a field 
which drives the system from an initial state at t = 0, 

\<p(t = 0)) = \<po), (3) 

to a target subspace at time T representing the final time 
objective, such that a minimum (or maximum) expecta- 
tion value of the time-dependent operator P(t) is main- 
tained throughout the complete time interval [0, T\. The 
target subspace at time T is described by the projector 
D, e.g., D = \ipf)(ff\ for a single target state. 



In OCT, these requirements are formulated as a func- 
tional which depends on the system state and the con- 
trol, in such a way that an optimal field corresponds to 
an extremum of the functional. That functional can be 
expressed as a sum over terms related to the different 
conditions imposed on the system evolution. 

1. The functional 

The complete functional is obtained as a sum over 
functionals corresponding to the final time objective, to 
the state-dependent intermediate-time objective (or con- 
straint), and to the constraint over the field. 

The term corresponding to the objective at the final 
time T, the actual target, can be expressed as 

Jo[<pr,<&] = \ (v(T)\D\<p(T)), (4) 

where Ao is a real parameter, which can be negative or 
positive, depending on whether the functional is mini- 
mized or maximized during the optimization. [t^T,^] 
emphasizes the bilinearity of the functional with respect 
to the system state at time T. Other possibilities for 
expressing this term exist [f| |24j , but the resulting opti- 
mization algorithms are very similar. 

The state-dependent intermediate-time objective or 
constraint is represented by the functional, 

Jb[<p,<p'] = / gb[f,<p']dt, (5) 
Jo 

where gb is taken to be 

g b [ip,^]=\b(v(t)\P(t)\<p(t)). (6) 

Ab is a real parameter which can be positive or negative, 
as discussed later. More complicated dependences of gb 
on the operator P(t) and on the state ip are conceivable. 

To obtain a closed algorithm, the complete functional 
has to include a term depending on the field 0, Q , 

Ja[e}= f g a [e]dt. (7) 
Jo 

Generally, g a can be written as 

g a [e]=X a (t)[e(t)-e r (t)] 2 , (8) 

where e r denotes a reference field and e r = corresponds 
to the common choice of minimizing the field energy. 
J a [e] represents an intermediate-time objective, but one 
which does not depend on the state of the system. 
The complete functional is given by 

J[<PiV^ e ] = Mvt^t] + Ja[e] + Jb[<P, y t ]. (9) 

For simplicity, we omit the dependence of ip and e on 
time, except for the final time T . The optimization prob- 
lem is now equivalent to the minimization or maximiza- 
tion of this functional. For that purpose, the Krotov 
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method is employed. Since the Krotov method operates 
with real functions, a complete presentation of the equa- 
tions for this problem is somewhat cumbersome [f3]. An 
outline of the derivation is given in the Appendix and 
only the final result is presented below. 



2. The optimization algorithm 

A guess field is denoted by (t) and the correspond- 
ing state \ip^(t)} is given by the evolution Eq. (p} with 
the initial condition \ip{t = 0)) = \ifo). In the Krotov 
method a new field eW(i) which decreases (or increases) 
the functional value is obtained by the following equa- 
tions: A new "state" \\) is determined using the inho- 
mogeneous equation 

d \x(t)) 



(/f ,v, ^H[e<°>(t)]| X (t)> + \ b P(t)\<pW(t)), (10) 

with the "initial" condition 

| X (T)> = -A D|^°)(T)), (11) 

cf. Eqs. (|A14IA17|) of the Appendix. It corresponds 
to the common OCT result modified by the inhomoge- 
neous term XbP(i)\ip(°\t)) which arises from the state- 
dependent constraint. The state \x) is used to determine 
the new control field, 



.(i) 



(*) 



(<) 



1 



-Im 



{<X(*)IAI¥> (1) (*)>} . (12) 



nx a (t) 

cf. Eq. (|A22|) of the Appendix, where e r = e^ -* was cho- 
sen. This is an implicit equation since the state 
which depends on appears in the right-hand side of 
Eq. p2[) . The numerical discretization of this implicit 
equation has been widely discussed for the homogeneous 
case (see for example Ref . @ ) . The inhomogeneous term 
in Eq. (I10|) requires a modification of the time prop- 
agation method. A symmetrical propagation scheme is 
employed based on the diagonalization of the Hamilto- 
nian in the interleaved time grid points, tt + At/2. The 
inhomogeneous term is evaluated as 



Xh 



P(t l+1 )\^\t l+1 )) + P(U)\^(U)) 



(0)1 



(13) 



The iterative algorithm is constructed with e^ 1 ' as in- 
put to the next step of the iteration and the process is 
repeated until the required convergence is achieved. 



3. Monotonia convergence 

The monotonic convergence of the algorithm is ana- 
lyzed defining A as the difference between the functional 
values before and after one iteration, 

A = Wt(o) >e (o)] _ j^W^tW^W] 

= A x + / (A 2a (t) + A 2b (t)) dt. (14) 



The terms Aj are derived in the Appendix and can be 
evaluated using Eqs. This yields 



Ai = -A <C(T)|D|C(T)>, 
cf. Eq. (|A"24|) . with the definition 

IC(t)> = M 1) (*)>-b (0) (*)> 

Furthermore, 

A 2Q (<) = -g a [e^]+g a [e^] + 

(«<-> 



(15) 
(16) 



dg a 



(0) 



(i) 



(17) 



cf. Eq. (|A29[) . which yields for our choice e r = 



A 2a (t) = X a (t) ( e (i)( t )_ e (0)( t )) a , (18) 



and 



A 2b (t) = -A 6 (C(t)|P(t)|C(t)>, 



(19) 



cf. Eq. (TA291) . 

The algorithm converges monotonically to a minimum 
(maximum) of the functional if A > (A < 0) in each 
iteration step. A sufficient but not necessary condition 
consists in all Aj being larger (smaller) than zero. Let the 
operators D and P(t) be positive-semidefinite. Sufficient 
conditions are then given by 



A <0, A 6 <0, A o (t)>0. 



for minimization, and by 



A >0, A b >0, A Q (t)<0 



(20) 



(21) 



for maximization. 

This result leads to some curious consequences. For 
example, let the system be described by a discrete num- 
ber of levels and assume its Hilbert space can be split 
into two subspaces, the "allowed" subspace, described by 
the projector P a iiow, and the "forbidden" subspace, de- 
scribed by P forbid (P allow + P forbid = 1)- The objective 
of the optimization consists in some transition inside the 
allowed subspace, avoiding any population transfer to the 
forbidden one. In the case of minimization of the func- 
tional J, the latter requirement can be expressed by one 
of the two following choices for J , 

(a) P(t) = P allow , X b <0, 

(b) P(t) = 

'forbid i 

X b > . (22) 
In the case of maximization, the possibilities are 

(a) P(t) = Pallow, X b > 0, 

(b) P(t) = P form, X b <0. (23) 

In both cases, (a) and (b) have the same physical mean- 
ing, remaining in the allowed subspace, or equivalently, 
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avoiding the forbidden subspace. The choice (6) is more 
appealing in principle, since the inhomogeneous term of 
Eq. (fTOf would decrease and eventually become negli- 
gible when approaching an optimal solution. However, 
only (a) fulfills the sufficient conditions for monotonic 
convergence. 

A note of caution must be made at this point. Eqs. 
(|2T))) and (|2"Tj) are sufficient but not necessary conditions. 
Monotonic convergence can therefore be found for values 
of A not fulfilling Eqs. flU and dHJ. This can happen 
if the values of Aj compensate each other to give a con- 
vergent total A. In addition, the analysis assumes an 
exact solution of the control equations. A limited accu- 
racy of the numerical implementation of the algorithm 
and a poor accuracy of the propagation method can lead 
to the breakdown of the monotonic convergence [26j . 



B. Optimization of a unitary transformation 

The objective consists in implementing a given unitary 
transformation O, up to a global phase, in a given sub- 
space Hp of dimension N r described by the projector R, 



(24) 



To this end, the parameter t is defined, 



Tr^U^O^R} = J2^f^n(T)), (25) 



71=1 



where 



0|n) = \ipfn), 
U(t,0;e)|n) = \<p n (t)) 



(26) 



The modulus of r is equal to N r when the target unitary 
transformation is implemented in the subspace 7i R by the 
field e @]. 

The optimization problem is again formulated as a 
functional minimization (maximization). The final time 
term is now defined by 



MWTn,^ Tn }} = A |r| 2 (27) 
= (<Pfn\<Pn(T)) E {<Pn>(T)\<Pfn>} , 



n=l 



where {ip n , <^t} denote the set of states \<Pn) (n = 
l,...,N r ). Other choices of Jo are possible p]. The 
intermediate-time state-dependent term takes the form, 

g b [W n ,^ n }} = A fc Tr{u(i,0; £ )tp(i)U(i,0;e)R} 

= \ b ^2(ip n (t)\P(t)\^ n (t)). (28) 



The constraint over the field is taken to be the same as 
in the state-to-state case, cf. Eq. (jHJ. 

The equations defining the algorithm are obtained us- 
ing the Krotov method as outlined in the Appendix. 
They read as follows: N r "states" \x) are given by the 
inhomogeneous evolution equation, 

d 'Xn(*)> = -^H[e<°>(i)]|xn(t)> + A 6 P|^(*)), (29) 



dt ,/v " y "" h 

with the "initial" condition 



\Xn{T)) = -A«> T> {6*0(2; 0;e<°>)ft}|p /n ) 

= -Ao ( £ (<Pf n >\<pW(T))) \Vfn), (30) 



n = 1, . . . , N r . They are used to determine the field 
by means of 

e«(t) = (31) 
e ( °H*)-^Im{E^WlAl^(*))}. 

The discussion of the monotony of convergence is 
equivalent to the state-to-state case @. That is, the suf- 
ficient conditions for minimization, Eq. (|20p . and maxi- 
mization, Eq. (|2ip are also valid for the optimization of 
a unitary transformation. 



III. ILLUSTRATION OF THE ALGORITHM 
ENFORCING A STATE-DEPENDENT 
CONSTRAINT 

In order to illustrate the algorithm outlined in the 
previous section, a simplified model of the vibrations 
in a Rb2 molecule is employed where only three elec- 
tronic states are considered, cf. fig. [TJ Of each elec- 
tronic state, 11 vibrational levels were chosen, specifi- 
cally v — 0, . . . , 10 from the A^E+ electronic ground 
state, v' — 5, . . . , 15 from the X E+ excited state, and 
v" = 2, . . . , 12 from the 1 II S excited state. A laser field, 
e(t), couples the X X T,+ levels to the 1 T,+ levels and the 
1 E+ to the 1 II g levels. In this model the transitions v = 
to v' = 10 and v' = 10 to v" = 6 have similar frequencies 
and Franck-Condon factors (FCF), cJi_>io = 0.0507 a. u. 
vs. tt>io^6 = 0.0506 a. u., and modulus of FCF 0.17 vs. 
0.23. The Hamiltonian describing our model is then given 
by 



H 



E H, ® \e l )(e l 



(32) 



n=l 



Ae(*) (|ei)(e 2 | + |e 2 )(ei| + |e 2 )(e 3 | + |e 3 )(e 2 |) , 

where a dipole moment operator fi independent of the in- 
ternuclear distance R is assumed. The vibrational Hamil- 
tonian is denoted by and the electronic state associ- 
ated to X 1 ^, and l H. g by \ei), (i = 1,2,3), re- 
spectively. The vibrational level energies and FCF were 
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FIG. 1: Schematic representation of the Rb2 model used in 
the calculations. The dotted lines indicate the position of the 
vibrational manifolds considered in each electronic state. 

obtained by diagonalization of the vibrational Hamilto- 
nians employing the potential energy curves of Ref . [2^ | . 

For the optimization examples described below, a guess 
field of the form 

e g (t) = e s{t) cos[fi(*-T)/2] (33) 

is employed. T corresponds to the target time for the 
optimization, set to T = 8ps, and s(t) — exp[— 32(t/T — 
1/2) 2 ] is taken to be a Gaussian shape function. The 
central frequency of the guess field is chosen to be f2 = 
u) v= o— >u'=io and eo = 10 -4 a.u.. The optimization param- 
eters related to the final time objective, Jo and to the field 
constraint, J a , are set to A = —1 and A a (t) = 100/s(i) 
in all calculations. 

The state-dependent constraint in the optimization 
forces the population to remain in the subspace of the 
two lower electronic states. This is formulated by iden- 
tifying P(i) = Fallow as the projector onto the X^Y^ 

and 1 S+ levels, P forbid corresponds thus to the projec- 
tor onto the 1 Ii g levels. This choice of the allowed sub- 
space is motivated as follows. The pulse duration is much 
shorter than the spontaneous emission lifetimes. On the 
timescale of the pulse, losses in a molecular system are 
due to processes such as predissociation, auto- or multi- 
photon ionization. Unlike spontaneous emission, these 
processes are relevant only in certain excited electronic 
states. Electronic states which are not affected by loss 
can therefore be included in the allowed subspace. 



A. State-to-state optimization under 
state-dependent constraints 

The objective for state-to-state optimization is cho- 
sen to transfer population initially in level v — of the 
electronic ground state to level v = 1 at time T, using 
Raman-like transitions via levels v' in the excited 
state. Optimizations with and without state-dependent 
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FIG. 2: (Color online) J n0 rm (violet solid line) and Ip (green 
dashed line) as a function of the number of iterations for \t = 
corresponding to optimization without state-dependent con- 
straint (upper panel) and XbT — — 32 (lower panel). Note the 
different scale of the z-axes in the two panels. 

constraint are compared. In both cases, the optimal field 
achieves a population transfer larger than 99.9% at the 
final time T. However, without the state-dependent con- 
straint the optimal field does not "know" that it is not 
supposed to populate levels in the upper electronic 1 IL g 
state at intermediate times. Due to similar transition fre- 
quencies and Franck-Condon factors, population transfer 
into the forbidden subspace therefore occurs. 

In order to compare the performance of optimization 
with and without state-dependent constraint quantita- 
tively, two measures are defined, 

Jnorm — T j > ™ j (34) 
A + Ab-l 

and 

Unlike the original functional Eq. ([9]), the normalized 
functional J nor m has an optimal value equal to one inde- 
pendent of the choice of Ao and Ah, while Ip corresponds 
to the average value of population in the allowed sub- 
space. Figure [5] shows the values of J norm and Ip as 
the iterative optimization proceeds. Optimization with 
state-dependent constraint requires a larger number of 
iterations to reach the optimal value of the total objec- 
tive J. However, this price is paid off since in this case 
the solution keeps indeed nearly all of the population in 
the allowed subspace at any time. 

Figures [3] and [4] demonstrate the behavior of the pop- 
ulations under the dynamics generated by the optimized 
pulses obtained with and without the state-dependent 
constraint. The amount of population in the upper elec- 
tronic state is largely reduced for a field obtained with 
state-dependent constraint as compared to that resulting 
from unconstrained optimization, cf. fig. [3J Two differ- 
ent transfer mechanisms are found: The pulse obtained 
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2 3 4 5 
time t [ ps ] 

FIG. 3: (Color online) Population in the forbidden subspace 
as a function of time for the system driven by the optimal 
field obtained with \ b = after 17 iterations (blue dashed 
line) and with XbT — —32 after 500 iterations (red solid line). 
The inset shows the population in the central time interval in 
a semi-logarithmic plot. 




FIG. 4: (Color online) Population in levels v = 0, v = 1, 
v' — 10 (in the allowed subspace) and v" = 6 (in the forbidden 
subspace) as a function of time for the system driven by the 
optimal field obtained with Ab = after 17 iterations (blue 
dashed line) and with XtT — —32 after 500 iterations (red 
solid line). Also shown is the population for the evolution 
with the guess field (black dotted line). Note the different 
scale of the y-axis for v" — 6. 



with the state-dcpcndcnt constraint transfers population 
to v = 1 in a ladder-like process which is driven by short 
subpulses. In between the subpulses, the amount of pop- 
ulation in the X E+ excited state (level v' = 10) is small, 
cf. fig. [U Further excitation to the forbidden 1 Ii g levels 
becomes thus unlikely. On the other hand, the pulse ob- 
tained without the state-dependent constraint transfers a 
large amount of population to the intermediate 1 S+ elec- 
tronic state (mainly to level v' — 10) which is later to be 
brought back to the ground state to level v = 1. The large 
amount of population which resides in the 1 £+ electronic 
state while the field is on allows for transient transfer to 
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FIG. 5: (Color online) Spectral amplitude as a function of the 
frequency u> for of the optimized field obtained with At = af- 
ter 17 iterations (blue dashed line) and with XtT = —32 after 
500 iterations (red solid line). The dashed vertical lines indi- 
cate the transition frequencies ui v=0 _ v t= 10 and u> v= i^ v '=io- 
The transition frequencies between the levels of the X E| 
electronic ground state and the 1 E„ intermediate state are 
represented by green circles, and the transition frequencies 
between the intermediate state and the 1 II 9 upper state 
by violet squares. 



the upper electronic state at intermediate times. 

Both transfer mechanisms share the common feature of 
a process driven mainly by one-photon transitions. Pop- 
ulation is transferred to the intermediate electronic state 
by a one-photon absorption. This population is later 
sent back to the electronic ground state by a one-photon 
emission or further excited to the upper electronic state 
by another one-photon absorption. Large spectral ampli- 
tudes of the optimal fields at frequencies corresponding 
to the main transition frequencies of the system reflect 
this finding, cf. fig. [5] Analysis of the two-photon spec- 
trum, i.e. of the Fourier transform of e(t)e(t), confirms 
that the amount of processes involving two or more pho- 
tons is very small. Figure [5] furthermore illustrates an 
important difference between the results of optimization 
with and without state-dependent constraint: The addi- 
tional requirement implies a more complex optimal solu- 
tion which is reflected both in a broader spectrum and 
in a more intricate dependence of the spectrum on fre- 
quency. 



B. Optimization of a unitary transformation under 
state-dependent constraints 

The implementation of a Fourier transform [30l | in lev- 
els v = 0, 1, 2, 3 of the X x £+ electronic ground state is 
chosen as example objective for the optimization of a uni- 
tary transformation. The state-dependent constraint is 
taken to be identical to the optimization of state-to-state 
transfer, i.e. the upper electronic state represents 
the forbidden subspace. Figure [S] shows the population 
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FIG. 6: (Color online) Population in the forbidden subspace 
as a function of time for the system driven by the optimal 
field obtained with \ b = after 50 iterations (blue dashed 
line) and with AbT — —8 after 500 iterations (red solid line) 
when the system is initially in level v = 0. 



in the levels of the forbidden subspace for the system 
driven by an optimized field obtained with and without 
state-dependent constraint. In both cases \t\ > 3.999; 
since the target value is equal to 4, the number of lev- 
els in which the unitary transformation is implemented, 
this corresponds to an error of less than 10 -3 . The main 
results are similar to those of optimizing a state-to-state 
transition: A larger number of iterations is needed to 
obtain similar performance with respect to the final time 
objective Jo. In addition the solution becomes more com- 
plex for the optimization with the state-dependent con- 
straint. The two tasks of a state-to-state transition and of 
a unitary transformation are run with the same number 
of iterations for optimization with the state-dependent 
constraint. As expected the state-to-state transition con- 
verges faster, cf. figs. [3] and [5] The effort for the opti- 
mization of a unitary transformation is approximately 
equivalent to N simultaneous state-to-state transitions 
(here N = 4), i.e. it corresponds to a more difficult opti- 
mization problem Q. 



C. Robustness with respect to decay 

A loss mechanism in the 1 II ff electronic state is mod- 
eled by adding an imaginary term — iT/2 to the vibra- 
tional energies, where T — 1/t£ denotes the decay rate 
and T£ the lifetime. Physically, such a decay is caused 
by processes such as predissociation or auto-ionization. 
When decay is included, the system dynamics generated 
by the optimized field is perturbed depending on the de- 
cay rate and on the amount of population in the lossy 
upper electronic state. Figure [7^ shows the population 
in the target level for the state-to-state transition with 
and without state-dependent constraint with optimiza- 
tion parameters as described above (solid lines). More- 
over, the total amount of population remaining in the 



FIG. 7: (Color online) Comparison of final results at time 
t = T as a function of the lifetime of the upper excited 
state, (a) State-to-state transition: Population in level v = 1, 
Pv=i(T), (solid lines) and total population remaining in the 
system, P S (T), (dotted lines) for the optimized field obtained 
with (red) and without the state-dependent constraint (blue) . 
(b) Unitary transformation: Normalized final time objec- 
tive, \r\/N r for a field obtained by optimization of the uni- 
tary transformation with (red solid line) and without state- 
dependent constraint (blue solid line). The vertical dotted 
line indicates the overall pulse duration. 

system at time T, 

P S {T) = &{T)\{P allow + P fo rbid)W{T)) , (36) 

is depicted by dotted lines in fig. [7^- The smaller popu- 
lation transfer to the upper state in case of optimization 
with the state-dependent constraint results in a larger 
robustness of the solution as the decay rate is increased. 
For a lifetime on the order of the pulse duration, the final 
time objective is only very slightly perturbed. In contrast 
the example without state-dependent constraint shows 
already a significant loss in the objective. This effect 
becomes even more evident when the lifetime is smaller 
than the pulse duration. For optimization without state- 
dependent constraint the transfer efficiency is reduced by 
50% for lifetimes on the order of 100 fs. For short life- 
times both optimization methods fail, nevertheless the al- 
gorithm including the state-dependent constraint is still 
superior. The decrease in efficiency of the transfer to the 
V = 1 target level (solid lines) is associated to the loss of 
population from the system (dotted lines) as the lifetime 
decreases. 

An interesting effect is found as the lifetime becomes 
very small: the final time objective actually improves. 
This is a manifestation of the quantum Zeno effect which 
states that a continuously observed quantum system 
never decays (3l| . In our example, the usual interpre- 
tation is inverted: The decay process can be associated 
with a weak measurement monitoring the population in 
the third electronic state, and the optimized field corre- 
sponds to the decay in the usual picture. As the decay 
rate increases, respectively the lifetime decreases, moni- 
toring of the forbidden subspace becomes continuous and 
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the pulse cannot populate the lossy state anymore. The 
possibility to loose population becomes thus smaller (cf. 
dottes lines in fig. [7K) which entails also better results for 
the final time objective (cf. solid lines in fig. [7^). How- 
ever, by comparing the solid and dotted lines in fig. [7^,, it 
is obvious that the loss of coherence in the control process 
is larger than the loss of population alone. 

Loss or decoherence is the main obstacle for generat- 
ing controlled unitary transformations. Figure [TJd shows 
how the objective deteriorates due to loss for the uni- 
tary transformation. The opimized field obtained with 
the state-dependent constraint is able to maintain a high 
fidelity despite the loss term. The performance of the 
optimized field for the unitary transformation is clearly 
worse than for the state-to-state transition. This is due 
to more population transfer to the upper electronic state 
in case of the unitary transformation obtained after the 
same number of iterations. A result comparable to that 
of optimizing the state-to-state transition with the state- 
dependent constraint can be achieved by increasing the 
number of iterations. 



IV. DISCUSSION 

The present work is related to a number of previ- 
ous OCT studies using state-dependent constraints, or 
equivalently, a time dependent target. The general for- 
mulation of our complete functional for optimization of 
a state-to-state transition employing a state-dependent 
constraint, Eq. ((9|), is related to the functional of Ref. 
[H| identifying D = X and P(£L= Y (t)- Similarly, 
the specific application of Refs. [2^, [2lJ corresponds to 
J = and P(t) = /(t)6; and the application of Ref. [HI 
to J = 0, P(t) = 6i(i) + 2TS(t - T)6 2 , and D = 6 2 . 
The approaches of Refs. [HI, Hy, [HJ consist in imposing 
a specific desired dynamics onto the system, for exam- 
ple, the populations following a given time profile, or the 
maximization of the expectation value of a given opera- 
tor. 

In the studies of Refs. [H, [H, [H, [27| , the optimization 
algorithm was obtained using the variational method. It 
is well-known that the variational approach is compatible 
with a large number of implementations of the control al- 
gorithm, not all of them showing monotonic convergence. 
Ref. [24| gives a detailed analysis of this family. Our re- 
sults which were obtained with the Krotov method, cor- 
respond to the case rjk = and = 1 of Ref. |24] . The 
Krotov method comes with the advantage of allowing a 
straightforward discussion of the convergence of the algo- 
rithm in terms of the sign of the optimization parameters 
Xi. As discussed in Section [Til this sign of the optimiza- 
tion parameter turns out to be crucial for the choice of 
projector P(t). 

As mentioned in the introduction, a different approach 
to avoid population leakage to undesired states is using 
local control theory 0, [2CJ, [U, [22j] ■ The example chosen 
in Ref. [22j] differs from ours: In Ref. (22j, the allowed 



subspace was restricted to the register levels of a unitary 
transformation. This choice forces the dynamics of the 
system to rely solely on two-photon processes which as- 
sure that the population remains in the register. A simi- 
lar task in our formulation would correspond to choosing 
the operator P allow as the projector onto the X 1 E+ levels 
only. However, the state-dependent constraint is formu- 
lated such as to maximize the time-averaged population 
in the allowed subspace. The family of solutions found 
by the algorithm consists in transferring some amount of 
population to the intermediate electronic state, but only 
for a very short time. This is reflected in the inset of fig. [3] 
in a sequence of "spikes" in the population of the forbid- 
den subspace as a function of time. Since the constraint 
of not populating any electronically excited state is much 
more restrictive than the time-averaged formulation, the 
number of optimization steps needed for reaching a spec- 
ified efficiency and hence the complexity of the solution 
increase largely compared to the present results. 



V. CONCLUSIONS 

Steering a system to a desired objective has to al- 
ways be balanced by the damage induced by the steer- 
ing process. In the present study optimal control theory 
was adopted to include a positive objective to be maxi- 
mized and a negative constraint to be avoided. To ensure 
monotonic convergence the additional constraint which 
depends on the state of the system was incorporated in 
the Krotov method. 

It turned out that the state-dependent constraint needs 
to be formulated in terms of maximizing population in 
the allowed subspace. While one could expect this to 
be equivalent to minimizing population in the forbidden 
subspace, the resulting algorithms and their convergence 
behavior differ markedly. 

The algorithm was applied to a simple model mimick- 
ing vibrational manifolds in three electronic states of an 
Rb 2 molecule. Population transfer from the vibrational 
level v = to v — 1 of the electronic ground state and 
the implementation of a Fourier transform in the levels 
v = 0,1,2,3 were chosen as optimization examples for 
a state-to-state transition and for a unitary transforma- 
tion. In both cases, the optimized fields induce Raman- 
like transitions via an electronically excited state. It was 
shown that optimization including the state-dependent 
constraint indeed avoids population of a higher electron- 
ically excited state. This state can either correspond to 
a loss channel itself or represent a resonant intermediate 
state in an unwanted multi-photon process. 

A similar task of promoting one objective while avoid- 
ing damage has been developed using local control theory 
20]. It is important to note that the optimal solutions 
are quite different. The local control scheme tries to min- 
imize the instantaneous population. Therefore it resorts 
to off-resonant two-photon transitions where only a tran- 
sient population exists. The OCT scheme minimizes the 
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integrated population in the forbidden subspace. As a 
result the solution can contain abrupt spikes of popula- 
tion. Since these spikes have no time to go anywhere the 
system remains protected against damage. 

The success of the scheme was demonstrated in the 
ability to cope with a real decay channel in the forbidden 
subspace. The decay term causes loss of population and 
loss of phase or decoherence. The state-to-state objec- 
tive can be thought to work without long term coherence. 
Brief periods of coherence are sufficient to generate the 
transitions. On the contrary the unitary operator objec- 
tive has to maintain phase coherence for the total period. 
The ability of the algorithm to find solutions that can 
cope with this scenario is encouraging. The robustness of 
the coherent control solution could result from a scheme 
based on a large number of interfering pathways. In this 
case loss of a few pathways will only slightly hinder the 
final objective. 
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which correspond to the total time derivative of tp(t) in 
Eq. ([1]). As mentioned above, the rigorous utilization of 
the Krotov method is somewhat cumbersome and little 
instructive. For simplicity the case of a functional de- 
pending on two real functions denoted by p and ip' and 
a real control e will be presented. This simple case can 
be connected to the original problem by the following 
relations, 



(A2) 



for ip = <p,x,f- Our derivation follows closely Ref. [||, 
but including the state-dependent constraint, i.e. the de- 
pendency on gb, cf. Eq. [S] The final equations for the 
original problem, Eqs. (| lOti 1 2[) . are obtained following 
the steps presented in this Appendix. Moreover, our out- 
line can be employed to derive new optimization schemes 
based on different choices of Jo, g a or gb- 



The scalar function $ and the functional L 



APPENDIX A: REVIEW OF THE KROTOV 
METHOD FOR OPTIMAL CONTROL THEORY 

In the following the system Hamiltonian H is assumed 
to be Hermitian. The derivation can easily be general- 
ized to non-Hermitian Hamiltonians [24j. Let's define the 
states /}, 



I/) 
(/I 



(<p\\&[e], 



(Al) 



All functions considered in the derivation, ip, e, etc., 
depend on t but for simplicity this dependence will only 
be made explicit for the initial and the final time, e.g. ipo, 
Pt- The Krotov method is based on the introduction of 
the arbitrary scalar function (p, ip') , the functions 



G(p T , <4, $t) = M<PT, <4) + *( T > <PT, <Pt) > (A3) 



and 



R(p, p\e, $) = -g a (e) - g b (p, p^) + 



dip 



f + r 



~dt 



(A4) 



r 



Note that 



R 



{9a+9b) + 



d<f> 

~dt 



where it was used that / (f') is the total time derivative 

of p (<£+). 

(A5) A new functional can be defined 
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'\<p,<p\e,$\ = G(<pt,¥t,®t) - $(0,^o,4) - / R(y,¥>\e,$)dt, 

Jo 



(A6) 



with the interesting property 

L[ip,ip\e,$\ = J[f,ip\e] 



(A7) 



The Krotov method takes advantage of this property and 
the freedom in the choice of $ to find an iterative algo- 
rithm that minimizes (maximizes) the original functional 



2. The iterative algorithm 



Let's start with a given field and the corresponding 
functions ip^ and (p^ ) . The values of the functional are 



denoted by J (0) and L (0) . The functional J (L) should 
be minimized (maximization will be discussed later on). 
The objective is then to determine a new control e^', 
given the functions ip^ and ip'si 1 ) ^ for which 



J(°) = L(°) > jW = L« 



(A8) 



For the time being, the control e and the functions ip, 
if' are treated as "independent variables" . The Krotov 
method accomplishes the optimization in two steps. 

1. A function $ is determined such that L[ip, f>\ e, <fr] 
has a maximum for <pt°> y ip^°> regardless of the 
control e. The maximum condition implies second 
order (functional) derivatives of L with respect to <p 
and <p' . Note that these derivatives should be eval- 
uated "at" any field e. However, the much simpler 
evaluation at e^ ^ turns out to be sufficient 

2. e^ 1 ' is determined such that the functional is min- 
imized with respect to all possible controls e. This 
implies a second order (functional) derivative of L 
with respect to e. Since $ is determined by the 
conditions of step 1, this can be done without con- 
sidering the effect of the change in the functional 
due to change in if and <pS induced by the new con- 
trol. 

Finally, the necessary relation between the field and the 
states through the evolution equation is imposed, result- 
ing in 



(A9) 



The conditions derived from the second order (func- 
tional) derivatives are rather complicated p|. Fortu- 
nately, for the problems under consideration, the min- 
imum and maximum conditions can be relaxed to ex- 
tremum conditions on the functional, for which only first 
order (functional) derivatives are needed. Moreover, the 



extremum conditions are common for minimization and 
maximization of the functional, therefore both cases are 
treated together below. Due to the relaxation of the min- 
imum and maximum to extremum conditions, an addi- 
tional step is required: The monotonic convergence to 
the target value of the final algorithm must be checked. 



3. First step: Determining $ up to first order 

The evaluation of an expression [. ..] at e^ ), <p(°> and 
i s denoted by [. . . ]( ). The extremum in L for ip(°> 
and ^t(o) corresponds to the following conditions on R: 



8R 



dip 



(0) 



OR 



(A10) 



(o) 



For example, using Eq. (|A4[) , the second derivative reads 



dR 



dip's 



'dip* 



9$ 

dip 



df_ + dp_ 

dip's dp^> 



"9<I>~ 


d 


"9<I>" 


d^ _ 


+ di 


_dip\ 



where we have used that 



(All) 



"99$" 


f + P 


d 9$" 


9 


' 9$ " 


d 


' 9$ " 


dip dip's _ 


dip's dip's _ 


+ dt 


_dipi_ 


~ It 


dip's _ 



New functions \ an d x^ arc defined, 
9$ 



(A12) 



x(t) 



dip's 



(0) 



9$ 



dip 



(A13) 



(0) 



The second condition in Eq. (|A10[) becomes then 



dx = _ j 
dt X 



df 



dip's 



(o) 



9/t 



dip's 



X 



(o) 



dgb 



dip's 



(A14) 



(0) 



and a similar result can be found for x ■ Moreover, the 
extremum condition on L also implies, 



dG 



dfr 



0. 



(0) 



dG 

dlf'rp 



0. 



(0) 



For example the second derivative gives 
dG dJ 9$ 



dip\ 



Using the previous definition of x, x(T) can be identified, 



dipt 



(A15) 



(A16) 



X{T) = XT 



dJ 

diptp 



(A17) 



J (0) 



11 



Therefore the extremum conditions allow to determine 
the function \ ( an d X s ) using the evolution Eq. (|A14|) 
with the "initial" condition at time t = T, Eq. (|A17|1 . 
With the knowledge of x ar >d X > the function $ can be 
constructed up to first order, 



(A18) 



The function <f>i of the original problem is obtained using 

Eq. dm, 

$i[t t <p,<pi] = ( X (t)\<p(t)) + (<p(t)\x(t)) 

= 2Re{( X (t)W))} . (A19) 



4. Second step: Determining e' 1 ^ 

The extremum condition on L with respect to the con- 
trol, evaluated in eW leads to 



dR 



de 



= 0. 



(A20) 



(i) 



where [-..jm denotes the evaluation of the expression 
[. . . ] at e' 1 ), <p^ and i^^ 1 ). Using the definition of R, 



the derivative gives 



dR 
de 



dg a 
de 



d<P 

dip 



df 
de 



a/t 

de 



5$ 



dip* 



(A21) 



where it was used that <I> does not explicitly depend on 
e. Using $ 1; the previous equation results in, 



dg a 
de 



(i) 



df 
de 



(1) 



a/t 

de 



X = 0. (A22) 



(i) 



Since <7 a and the functions / and f\ defined in Eq. (|Al|l . 
depend on the control, this equation allows to determine 



5. Conditions for monotonic convergence 

To check whether the iterative method shows mono- 
tonic convergence, we define Ai and A 2 {t), 



Lfo<V t( V 0) ,*i 0) ] " Lfc> (1 V (1) ,e« ^ = A 1+ / A a (t)dt 



(A23) 



where 



To simplify this expression, it is assumed that 



= j8 S> -*®-# ) +*®, (A24) 



(!) „t(l) 



a/t 

dip 



= 0, 



a/ 
a<^t 



= 0, 



(A26) 



and 



A 2 (t) = R{ip {1) , </? t(1) , e (1) , - i?<y o) , </? t(0) , e (0) , $J 0) ) . which is true for the original problem, and the function 

(A25) R is split into R = R a + R b , 



R a = -g a - X* 



a/ 

dip 



<P + X f / - t P t 



(o) 



a/t 



(0) 



-gb 



dgb 



dip 



ip + cp 



(o) 



(A27) 
(A28) 



(0) 



Therefore A 2 (t) = A 2a (t) + A 2b (t) with 
A 2a (t) = [i?a] ( i) - [Ra](o) 

A 2 &(t) = [i?b] (1) - [-R&](o) 



(A29) 



If the objective is to minimize the functional J (j( Q ) > 
J^), the sufficient but not necessary conditions for 
monotonic convergence are given by Aj > 0, (j = 1, 2, 3). 
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Analogously, in order to maximize the functional it is suf- 
ficient that Aj < 0, (j = 1,2,3). Whether the Aj are 
positive or negative is determined by the particular choice 
of Jo, <7a, and 9b- In Sec. [IT] the values of A arc analyzed 
for the cases under study. 

6. The Krotov method for unitary transformations 

The previous approach can easily be generalized to 
the case of unitary transformations. The functional J 
depends then on 2N r real functions denoted by the set 
{<Pn,<Pn}- As a consequence, the scalar functions G 
and R will also depend on all of them. Moreover, the new 



dependence must be taken into account in the derivatives 
of the previous equations by the substitution 

d_ V — 

w - ?£■ 

The remaining procedure is analogous to the state-to- 
state case, and the relations (|A2|) can be used to ob- 
tain the equations for the optimization algorithm given 
in Sec. HTbI 
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